Cosmological hydrogen recombination: The effect of extremely high-n states 
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Calculations of cosmological hydrogen recombination are vital for the extraction of cosmological 
parameters from cosmic microwave background (CMB) observations, and for imposing constraints 
to inflation and reionization. The Planck mission and future experiments will make high precision 
measurements of CMB anisotropies at angular scales as small as t ~ 2500, necessitating a calcula- 
tion of recombination with fractional accuracy of ~ 10~ 3 . Recent work on recombination includes 
two-photon transitions from high excitation states and many radiative transfer effects. Modern 
recombination calculations separately follow angular momentum sublevels of the hydrogen atom to 
accurately treat nonequilibrium effects at late times [z < 900). The inclusion of extremely high-n 
(n > 100) states of hydrogen is then computationally challenging, preventing until now a determi- 
nation of the maximum n needed to predict CMB anisotropy spectra with sufficient accuracy for 
Planck. Here, results from a new multi-level-atom code (RecSparse) are presented. For the first 
time, 'forbidden' quadrupole transitions of hydrogen are included, but shown to be negligible. Rec- 
Sparse is designed to quickly calculate recombination histories including extremely high-n states 
in hydrogen. Histories for a sequence of values as high as n max = 250 are computed, keeping track 
of all angular momentum sublevels and energy shells of the hydrogen atom separately. Use of an 
insufficiently high n max value (e.g., n max = 64) leads to errors (e.g., 1.8a for Planck) in the predicted 
CMB power spectrum. Extrapolating errors, the resulting CMB anisotropy spectra are converged 
to ~ 0.5ct at Fisher-matrix level for n max = 128, in the purely radiative case. 

PACS numbers: 98.70.Vc,32.70.Cs,32.80.Rm,98.80.-k 



I. INTRODUCTION 

Measurements of cosmic microwave background 
(CMB) temperature anisotropies by the Wilkinson Mi- 
crowave Anisotropy Probe (WMAP) have ushered in the 
era of precision cosmology, confirming that the Universe 
is spatially flat, with a matter budget dominated by dark 
matter and a baryonic mass fraction f2f,/i 2 [l| in agree- 
ment with the measured ratio of deuterium-hydrogen 
abundances (D/H) Q. WMAP measurements of large- 
scale CMB polarization also yield the optical depth t to 
the surface of last scattering (SLS), meaningfully con- 
straining cosmological reionization. Together with sur- 
veys of supernovae [!, H| , galaxies 0-0] , and galaxy clus- 
ters (H, WMAP measurements build the case that the 
Universe's expansion is accelerating, due to "dark en- 
ergy" or modifications of general relativity d, [l(J, and 
constrain other physical parameters (such as the sum of 
neutrino masses m Ui |llTjl3| and the effective number 
of massless neutrino species N v ) . 

CMB temperature observations (WMAP, 



BOOMERANG |14j, CBI [15| and ACBAR [16|) 



probe properties of the primordial density field, such as 
the amplitude A s , slope n s , and running a s of its power 
spectrum. These observations constrain deviations 
from the adiabatic, nearly scale free and Gaussian 
spectrum of perturbations predicted by the simplest 
models of inflation, but also offer controversial hints 
of deviations from these models (see Refs. 0, [I?} 
and references therein). Experimental upper limits to 
B-mode polarization anisotropies (e.g. DASI [l8[ and 
BICEP pJl) impose constraints to the energy density of 
relic primordial gravitational waves [l^, [2l[ . 



The Planck satellite, launched in May 2009, will obtain 
extremely precise measurements of the CMB tempera- 
ture anisotropy power spectrum (Cj T ) up to i ~ 2500 
and the E-mode polarization anisotropy power spectrum 
(Cf E ) up to t ~ 1500 [I!]. Robust measurements of 
the acoustic horizon and distance to the SLS will break 
degeneracies in dark energy surveys 0, I22T4241 ) . Polar- 
ization measurements will yield the optical depth r to 
the SLS [13], further constraining models of reioniza- 
tion and breaking the degeneracy between n s and r [22| . 
Cosmological parameters will be determined with much 
greater precision. More precise values of n s and a s will 
be obtained from CMB data alone, helping to robustly 
constrain inflationary models and alternatives to infla- 
tion [221. The advent of Planck, ongoing (SPT (2f| and 
ACT [20|) experiments at small scales, and a future space 
based polarization experiment like CMBPol [13, HH all 
require predictions of primary anisotropy multipole mo- 
ments Ci with 0(1O -3 ) accuracy. 

During atomic hydrogen (H) recombination, the 
Thomson scattering opacity drops, decoupling the 
baryon-photon plasma and freezing in acoustic oscilla- 
tions. The phases of acoustic modes are set by the peak 
location of the visibility function [29|, [3(| , while d amping 
scales [U Hl| and the amplitude of polarization [33[ H3| 
are set by its width. Small-scale CMB anisotropies are 
also smeared out by free electrons along the line of sight, 
suppressing power on small scales so that Ce — > Cge~ 2T , 
where r is the total optical depth of this w [35}. An 
accurate prediction of the time-dependent free-electron 
fraction x e (z) from cosmological recombination is thus 
essential to accurately predict CMB anisotropies. 

Recent work has highlighted corrections of 
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Ax e (z)/x e (z) > 0.1% to the standard recombination 
history computed by RecFast 36]. These corrections 
will propagate through to predictions of anisotropics, 
and neglecting them would lead to biases and errors 
in Planck measurements of cosmological parameters 
[13, EH- The use of the CMB as a probe of the first 
ionizing sources and of physics at energy scales greater 
than 10 16 GeV thus requires an accurate treatment of 
the ~ eV atomic physics of recombination [39| . 

Direct recombination to the hydrogen ground state is 
ineffective because of the high optical depth to photoion- 
ization [Z(| EH . Recombination proceeds indirectly, first 
through recombination to a n > 2 state of H, and then 
by cascades to the ground state. Because of the optical 
thickness of the Lyman-n (Lyn) lines, the resulting radi- 
ation may be immediately absorbed, exciting atoms into 
easily ionized states. 

There are two ways around this bottleneck [Z(], |4l| . In 
the first, the sequence of decays from excited H levels 
ends with a two-photon decay (usually 2s — > Is). The 
emitted photons may have a continuous range of ener- 
gies, allowing escape off resonance and a net recombina- 
tion. In the second, photons emitted in the np — > Is 
transition redshift off resonance due to cosmological ex- 
pansion, preventing re-excitation and yielding some net 
recombination. The dominant escape channel is from the 
2p — > Is Lyman-a line. These resonant transitions give 
off line radiation and distort the CMB 0, El]. 

Peebles, Sunyaev, Kurt, and Zel'dovich modeled re- 
combination assuming that all net recombination re- 
sulted from escaping the n = 2 bottleneck [Z(], |4l| . This 
three-level-atom (TLA) treatment included recombina- 
tions to excited states, under the assumption of equi- 
librium between energy levels n and angular momentum 
sublevels I for all n > 2 (note the use of I for atomic angu- 
lar momentum and £ for CMB multipole number). This 
sufficed until the multi-level-atom (MLA) model of Sea- 
ger et al. [3], which included hydrogen (H) and helium 
(He), separately evolved excited states assuming equi- 
librium between different I, accurately tracked the mat- 
ter/radiation temperatures Tm/Tr, [45], |46[, accounted 
for line emission using the Sobolev approximation [47j . 
and included H2 chemistry. This treatment underlies the 
RecFast module used by most CMB anisotropy codes, 
including those used for WMAP data analysis |3q ]. 

The higher precision of Planck requires new physical 
effects to be considered, among them two-photon tran- 
sitions from higher excited states in H and He |48l- l 52 | . 
other forbidden and semiforbidden transitions in He [53l — 
[55| , feedback from Lyn lines [5(| , and corrections to the 
Sobolev approximation due to a host of radiative trans- 
fer effects in H and He resonance lines [H, . Most 
recent work on recombination has focused in one way or 
another on the radiative transfer problem. Here we direct 
our attention to the populations of very high- n states. 

One important effect is the breakdown of statistical 
equilibrium between states with the same value of the 
principal number n but different angular momenta I. 



This effect is dramatic at late times. When I sublevels 
of a level n are resolved, increases in x e (z) of ~ 1% at 
late times result This changes predicted C/s at 

a statistically significant level for Planck. Highly excited 
states in hydrogen also change the recombination history 
at a level significant for Planck. While levels as high 
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as n = 300 were included in the treatment of Ref. [44 1 
underlying RecFast, I sublevels were not resolved. It 
is thus important to update cosmological recombination 
histories to include high-n states of H while resolving I 
sublevels, in order to predict the Cgs as well as CMB 
spectral distortions from recombination. 

Simultaneously including very high n and resolving the 

1 sublevels is computationally expensive, taking nearly a 
week on a standard workstation for n max = 100 |6l|, us- 
ing a conventional multilevel-atom recombination code. 
This becomes prohibitively expensive for higher values 
of n max , unless considerable resources are devoted to the 
problem. To date, this has prevented a determination of 
how x e (z) converges with n max and how high n max must 
be to predict CV's for Planck. The existence of electric 
dipole selection rules Al = ±1 means the relevant rate 
matrices are sparse, and we have used this fact to develop 
a fast code, RecSparse, to explore convergence with 
n m ax- While the computation time t comp for standard l- 
resolving recombination codes scales as t, 
with RecSparse the scaling is i comp oc n° 

2 < a < 3. With RecSparse, we can calculate recombi- 
nation histories for n max = 200 in 4 days on a standard 
work-station; this would likely take weeks using a con- 
ventional code. For the first time, we have calculated 
recombination histories for n max as high as 250 with I 
sublevels resolved. 

While previous computations have included some for- 
bidden transitions, none have included optically thick 
electric quadrupole (E2) transitions in atomic hydrogen. 
We include E2 transitions, and find that the resulting 
correction to CMB anisotropies is negligible. 

We find that the correction to CMB C/s due to ex- 
tremely excited levels is 0.5er or less if n max > 128, in the 
purely radiative case. This paper is not the final word 
on recombination; atomic collisions must be properly in- 
cluded and the effect of levels with n > n max must be in- 
cluded to conclusively demonstrate absolute convergence. 
The end goal of the present recombination research pro- 
gram is to include all important effects in a replacement 
for RecFast, as the interplay of different effects is sub- 
tle. 

In Sec. HH we review the formalism of the multilevel 
atom (MLA), and follow by explaining how we extend 
the MLA to include very high-n states (Sec. IIII[) and 
electric quadrupole transitions (Sec. HVp . State popula- 
tions, recombination histories, and effects on the CV's are 
presented in Sec. [V] We conclude in Sec. [VT] 

We use the same fiducial cosmology as in Ref. f62| : 
total matter density parameter Q m h 2 = 0.13, ft^h 2 — 
0.022, T C mb = 2.728 K, N v = 3.04, and helium mass 
fraction Y Rc = 0.24. 
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II. THE STANDARD MULTILEVEL ATOM 

We now review the elements of the standard multilevel- 
atom (MLA) treatment of cosmological recombination. 
For fundamental constants, we use NIST (National Insti- 
tute of Standards and Technology) CODATA (Commit- 
tee on Data for Science and Technology) values every- 
where 63]. Unless explicitly noted otherwise, we make 
the substitution m e — > fj, = m e m p / (m e + m p ) in all ex- 
pressions for the Bohr radius clq and the ground-state 
hydrogen ionization potential I R to correctly account for 
reduced-mass effects. 



where N (E, T R ) is the photon occupation number at 
photon energy E and radiation temperature T R . Here 
e is an infinitesimal line width and E n>nt is the energy of 
a photon produced in the transition [n',Z'] — > [n,l]. The 
simplest possible assumption for Af(E,T R ) is a black- 
body; we discuss further subtleties in Sec. Ill Bl 



N(E n 



,T R ) 



1 



,/(kT Ti ) 



(4) 



Here k is the usual Boltzmann constant. The (l +J\f^ n ,) 
term accounts for stimulated and spontaneous emission. 
The two-photon term is (ioH iHsij 



A. Basic framework 

CGS units are used except where explicitly noted oth- 
erwise. We follow the abundance x n ,i — ?7n,z/?7H, where 
?7h is the total number density of hydrogen nuclei and 
rj n> i is the density of hydrogen in a state with principal 
quantum number n and angular momentum I (we denote 
the state [n,Z]). We evolve these abundances including 
bound-bound and bound-free radiative, single photon, 
dipole transitions, as well as the 2s — > Is two-photon 
transition, which has rate A 2s ,i s = 8.2245809 s _1 [gj]. 
Focusing on the effect of single-photon dipole processes at 
high n max , we neglect higher n two-photon processes but 
note that their effects are large enough that they must be 
included in a final recombination code [48T - l5ll l53l | . Note 
that we also neglect collisional transitions. We comment 
on how this may change our conclusions in Sec. IV A II 

Bound-bound electric dipole processes are described by 
the equation [H H3 



(r M \ 

\ n,n 



•En' I' 



■pi', I 
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with 
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<>i(l+A/l<) ifn'>n, 
{ A%Pl/ n {gi/g v )N+ n if n'<n, 



(1) 



(2) 



where A' , is the downward Einstein rate coefficient for 

decays from [n',l'] to [n,l] and f„' n / is the probability 
that a photon emitted in the [n',l'] — > [n,l] line escapes 
the resonance without being reabsorbed. This probabil- 
ity is calculated in the Sobolev approximation, described 
in Sec. Ill Bl For lower / states easily described using 
the s,p,d,f... orbital notation, we will sometimes use 
the notation A®'^ = A nPt i s , Pi' n = P n p,is, and so on 
to simplify the discussion. The degeneracy of [n, 1} is 
gi = 2(21 + 1). We explicitly keep track of the angu- 
lar momentum quantum number I, as this will simplify 
discussion of our sparse-matrix technique in Sec. IIIIBI 

The photon occupation number blueward/redward of 
a line transition ([n',Z'] — > [n,l]) is denoted 

Kt=^(En^±e,T R ), (3) 



A, 



^2s->ls 



-x 2s 



2"; 



xi s e 



27 

■E 2Btls /(kT R ) 



(5) 



where E2 S ,i s = Ei,\ and the second term describes two- 
photon absorption with a rate coefficient obtained by re- 
quiring that forward/backward rates satisfy the principle 
of detailed balance. 

The bound-free term is E3, H [13 



bf 



= J[m 



x\a n i (E e ) S - x n jl (E e ,T T )] dE e , (6) 



with 



S(E e ,T M ,T R ) = [l+Af(E 7 ,T R )}P M (E e ,T M ) (7) 



and 



I(E e ,T R )=p nl (E e )Af(Ey,T R ) 



(8) 



This integral is over the total energy E e of a recombin- 
ing electron. The energy of a recombination photon is 
E 1 = E e — E n , where E n is the bound-state energy of the 
recombined electron. The recombination rate in cm 3 s _1 
of such an electron to the bound state [n, I] is a n i (E e ) 
and is discussed in Sec. IIII A 21 The ionization rate in 
s _1 is Pra (E e ), and easily shown by detailed balance con- 
siderations to be I5ll 



p nl (E e ) = a nl (E e ) 



2 7 /V/E^ 
h 3 gi 



(9) 



The free-electron abundance is x e = r) e /rjn, where r] e 
is the free-electron density. We restrict our attention to 
times after helium recombination, and so the free proton 
abundance x p — x e . The net bound- free rate [Eq. ©] in- 
cludes both spontaneous and stimulated recombination. 
The electron energy distribution is a Maxwellian with 
matter temperature Tm- 



P M (T M ,E e ) = 2. 



E R 



,-E e /{kT M ) 



7r(fcTi 



(10) 
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B. Radiative transfer and escape probabilities 

Numerically solving the radiative transfer problem is 
computationally intensive, but tremendous simplification 
can be achieved with the Sobolev escape probability for- 
malism, also known as the Sobolev approximation [47| . 
The Hubble flow can be used to define a lengthscale over 
which the bulk flow induces a velocity change equal to 
the thermal velocity: L = y/ '3kTM/ 'm a tom/ H (2r) , where 
H (Tr) is the value of the Hubble expansion parameter 
when the radiation has temperature Tr and m a tom is the 
mass of an atom 14411 . The conditions of the Sobolev ap- 
proximation are [44l . |47| . |62| : (i) L is much smaller than 
the typical length scales over which cosmological quanti- 
ties vary, (ii) L/c is much smaller than the typical time 
scales over which cosmological quantities vary, (iii) com- 
plete frequency distribution — the rest-frame frequency 
of an outgoing scattered photon v does not depend on 
the incoming frequency v' — and (iv) no other emission, 
absorption, or scattering processes occur in the vicinity of 
the line. Corrections to the Sobolev approximation result 
from diffusion around resonance lines [65l . |66| , atomic re- 
coil [H,[63], Thomson scattering near resonances [H,[6^], 
and overlap of the higher Ly series lines, leading to im- 
portant corrections to cosmological recombination calcu- 
lations. In this work, however, we work in the Sobolev 
approximation to focus on other physical effects. 

In the Sobolev approximation, the escape probabil- 
ity for photons produced in the downward transition 
[n', I'} -> [n, 1} is 



P 



1 



i.v 



where the Sobolev optical depth is given by 



i,i' 



*«B<n> '-■ V'// 



c m * i,i' / gi' 

^r, n> I X n,l ~ X n l V 



with transition frequency 

_ E nin ' _ In 



1 1 



(11) 



(12) 



(13) 



Correct expressions for n' < n are obtained by reversing 
arguments. During cosmological recombination, transi- 
tions between excited states are optically thin {P^ l n i > 

0.99972) and so we set P l / n , = 1 in our calculations 
for non-Lyman lines. 

Transitions in the Lyman (Ly) series (n' > n = 1, 
I' = 1, I = 0) are optically thick (r^, » 1) 



and so P^'^j — lAi' n '- Ly transitions cannot, how- 
ever, be ignored in the recombination calculation, as the 
rate at which atoms find their way to the ground state 
through the redshifting of resonance photons, P®'*, JL \ , 
is comparable to A2 S ->i s and other two-photon rates [51j , 
Strictly speaking, r°'^/ depends on and so one 



should solve for x n /,i and then iteratively improve the 
solution. The populations of the excited states, however, 
are very small and the maximum resulting correction to 
the optical depth is 2 x 1CT 12 (for ri =2,z= 1600) [H|. 
We thus drop the second term in Eq. (|12l) . simplifying 
our computation by working in the approximation where 
the Lyman-n (Lyn) line optical depth depends only on 
the ground-state population and not on the excited-state 
populations. 

Another aspect of the Lyman-series lines is feedback: 
a photon that escapes from the Lyn (np — > Is) line will 
redshift into the Ly(n — 1) line and be reabsorbed. Rec- 
Sparse has the ability to implement the resulting feed- 
back, using the iterative technique of Ref. TO]. This 
slows down the code by a factor of a few, however, and 
so to efficiently focus on the n max problem, we turned 
feedback off. For the high Lyman lines, feedback is al- 
most instantaneous: the Universe expands by a factor of 
A In a ~ 2n~ 3 during the time it takes to redshift from 
Lyn to Ly(n — 1). In the instantaneous-feedback limit, 
the Lyn lines do not lead to a net flux of H atoms to the 
ground state. To approximate this net effect we turned 
off Lyman transitions with n > 3; this leads to a smaller 
error than would result from leaving these transitions on 
but disabling feedback. Previous tests using the code 
of Ref. [(32| show resulting errors in the recombination 
history at the « 1% level; in any case, this should only 
weakly be related to the n max problem. All of the recom- 
bination histories and plots in this paper were produced 
by running RecSparse with both feedback and Lyman 
transitions from n > 3 disabled. 

Electrons, though nonrelativistic during recombina- 
tion, interact with photons through Thomson scattering. 
As a result, they do not follow the simple adiabatic scal- 
ing Tm oc a~ 2 , where a is the cosmological scale factor. 
The matter temperature is set using the asymptotic solu- 
tion of Ref. [5l| for z > 500, after which the relevant or- 
dinary differential equation (ODE) is solved numerically; 
this transition point occurs in the regime of mutual va- 
lidity for the numerical and asymptotic solutions. We ne- 
glect subdominant processes, such as free-free, line, pho- 
torecombination and collisional ionization cooling, as well 
as photoionization and collisional recombination heating 



as pr 
0. 



C. The steady-state approximation 

The wide range of disparate time scales in this prob- 
lem would naively necessitate a stiff differential equation 
solver. This computational expense can be avoided by 
repackaging Eqs. JT]), @,and ([S])-©. These equations 
may be rewritten for excited states as Qn, I] ^ [1,0]) 



•En,l 



n'V 



l.l' 



(14) 
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with 



T n% = 5 n'n' + 7nJ + J2 T n'L' ~ T n,n>i ( 15 ) 



n" ,1" 



where the integrated photoionization rate from [n, I] is 
2„l = J Pnl (E e ) I (-B e , T R ) dE e (16) 

and is defined in Eq. 

The downward flux to the ground state is 



lnl = Al%P?£ (1 + A/-+) Si tl + A 2s , ls £ 



(17) 



where the first term describes Lyn series transitions 
(stimulated and spontaneous) while the second accounts 
for the [2,0] — > [1,0] two-photon transition. Kronecker 

delta symbols (5^ n , and <5y) are employed throughout 
to enforce [n, 1} — [n', I'] and I = I'). 

The source term s n i includes flux from the ground state 
and direct recombination into the state [n, I] : 



s n ,l 



rjnxl J U n l (E E ) S (E e , T M , T R ) dE e 

°n,2 



xi s A 2Si i s e E2s 

|0,1 n 0.1 



x ls9l A^P^Af+S lA /2. 



(18) 



This can also be rewritten in matrix notation: dx/dt = 
—Tx + s, where T is the matrix of rates with components 
given by Eq. (fl~5|) . 

The left-hand side of Eq. (fH)) is associated with the 
recombination time scale, while both terms on the right- 
hand side are associated with much shorter atomic time 
scales. For example, the longest lifetimes in the re- 
combination problem are those of the 2s and 2p states 
(A-2 S ,is ~ 10 s and A2 Pt i s P2p,is ~ 1 s when Ly-a opti- 
cal depth peaks at r ~ 6 x 10 8 ), considerably shorter 
than the recombination time scale of t rec ~ 10 12 s. Thus 
we make a steady-state approximation, x n i = 0, which 
is formally valid because the reciprocal of the minimum 
eigenvalue of T peaks at 0.8 s, which is ~ 10~ 12 of the 
duration of recombination. Thus the excited-state abun- 
dances are given by 



(19) 



The rates in T and s depend on x e , x\ s , Tm, Tr, and 
J\f. The ground-state population is given by x\ s = 1 — 
x e — J2[ n i]jt[i o] x n,u but since excited-state populations 
are small (x n j < 10 -13 ), x\ s can be eliminated from 
Eq. (fT5|) using the approximation x\ s ~ 1 — x e . We can 
then solve for the evolution of x e , leaving out ineffective 
direct recombinations to the ground state: 



~ -±u = x ls A 2s ,ise- E ^/^ 

E[n,i]^[l,0] nnlX n ,l ~ T^hnPhn^lnXlsSl,! 



(20) 



The steady-state approximation thus allows us to con- 
vert a stiff system of ordinary differential equations into a 
large system of coupled linear algebraic equations, along 
with a single ordinary differential equation. 



III. RECOMBINATION WITH HIGH-N STATES 

The original "effective 3-level atom" (TLA) treatments 
of cosmological recombination in Refs. [40,|41| were built 
on the assumption that the primary bottlenecks to effec- 
tive recombination are the slow 2s — > Is transition rate 
and the reabsorption of 2p — > Is resonance photons by 
the optically thick plasma. Other crucial assumptions 
included radiative equilibrium between excited states, 

Xn = ^ 2e -(^-^)/( fer «)n 2 /4 if n > 2, (21) 

Xn = ^ y x n i 7 (22) 
l<n— 1 

and statistical equilibrium between angular momentum 
sublevels: 



Xn l 



(21 + 1) 



(23) 



Recombination to higher excited states was included 
through an effective "Case B" total recombination con- 
stant (T) (recombinations to the ground state are 
omitted) gO, H. 

As the radiation field cools and the baryon density falls 
at late times, the transitions coupling high-n to low-n be- 
come inefficient, as do those coupling different / sublevels 
with the same n. This leads to a breakdown of statistical 
equilibrium (note however that the steady-state approx- 
imation is still valid), and so Eqs. (|2"2")l and (|2"3")l cease to 
apply. In Ref. flij]. Eq. (|2"2"]l is relaxed while Eq. $Sfr 
is still imposed, and ~ 10% corrections to the TLA pre- 
diction for x e (z) result. At late times, nonequilibrium 
effects cause a net flux downward from states with quan- 
tum number n to the ground state, accelerating recom- 
bination. 

The inclusion of progressively more shells increases the 
number of downward cascade channels to the ground 
state for continuum electrons. Thus higher n max leads 
to faster recombination and lower x e (z). Reference [44| 
reports results for n max as high as 300. The Lyman 
(np — > Is) transitions from very high-n states overlap 
with the Lyman continuum, motivating Ref. [44j 's claim 
that there is no need to go past n = 300. The real ques- 
tion as to whether the different values of n are well de- 
fined, however, is whether the broadening of the state, 
H/t (where r is the lifetime) is larger than the splitting 
of adjacent energy levels, AE s» 2Inn~ 3 . The intrin- 
sic broadening for a typical level with l/n ~ 0(\) is 
H/t ~ a 3 / H n -5 (z3- Thus H/t « AE and so these 
extremely high-n energy levels are well defined; indeed, 
transitions between highly excited states in such nonequi- 
librium plasmas are seen in interstellar H 11 regions and 
are a useful diagnostic of physical conditions [721 ] . 
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For extremely large n, the above physical argument 
may break down because of additional broadening con- 
tributed by interactions with the radiation field and the 
plasma. For example, the broadening due to stimulated 
emission and absorption scales as ~ n~ 2 (the sponta- 
neous n~ 5 times the phase space density for photons 
in the An = ±1 transitions) and that due to electron- 



impact collisions scales as 



[731 ]: at sufficiently high 



n these will dominate over n -3 and the atomic energy 
levels will become blended. However, the orders of mag- 
nitude of the collisional coefficients [73[ suggest that this 
occurs at values of n larger than those considered in this 
paper. We have also verified that for conditions of inter- 
est for the recombining cosmological plasma, the plasma 
Debye length is greater than the average bound electron 
radius a^n 2 as long as n £ 10 5 . 



More recent work [60|,|6l| shows that additional ~ 1% 
corrections to x e {z) arise when Eq. (|23l) is not imposed 
and the populations of I sublevels are followed separately. 
Bottlenecks to decays from high I imposed by V = I ± 1 
slow down cascades to the ground state, and thus lead to 
slower recombination. In this case, the sidelength of T is 
N = O (n^J . Since the number of computational steps 
needed to invert a matrix is generically a N 3 process, the 
computational time needed for a single ODE time step 
in the recombination time will be proportional to 



i 6 

max • 



As noted in Ref. [6l| , a recombination calculation with 
n max = 100 already takes ~ 6 days on a standard work- 
station. It this thus difficult to explore how quickly x e (z) 
converges for progressively higher values of n max . Even 
between n max = 80 and n max = 100, ~ 1% changes are 
seen in the TT and EE multipole moments (C/s) of the 
CMB 1 61]. In spite of the computational challenge, it is 
thus crucial to push the calculation to sufficiently high 
"max that corrections to x e (z) from remaining n > n m ax 
are so small that they do not effect Cj T or Cf E at a level 
statistically significant compared to the predicted Planck 
sample variance (e.g., several parts in 10 4 for I > 1000) 
(69j . There are two challenges in treating such a big mul- 
tilevel atom. The first is the calculation of atomic transi- 
tion rates at extremely high n; this is tractable because of 
some convenient recursion relations. The second is simul- 
taneously evolving the populations of n max (n max + 1) /2 
states. We discuss these in turn below. 



1 In Ref. [74j . the results of Ref. [6l| are used to explore the ef- 
fect of progressively higher n max on CMB Cg's. In that work, 
It is noted that the fractional difference between the C/s for 
"max = 60 and ra max = 120 falls within a heuristic Planck per- 
formance benchmark. Higher values of ti max come even closer to 
the fiducial case of n max = 120, a fact used to argue that even 
"max = 60 recombination is adequate for Planck data analysis. 
From the Cauchy convergence criterion, however, we know that a 
meaningful convergence test requires a comparison between suc- 
cessive members in a sequence. Using the results of Ref. [6l| 
alone, the question of convergence with n max thus remains open. 



A. Rates 

Here we discuss the Einstein coefficients for dipolc 
bound-bound and bound-free transitions in atomic hy- 
drogen, which are used in our recombination computa- 
tion. We omit reduced-mass corrections to make a con- 
sistent comparison with Refs. [ZlllzB^Zll; but include 
them when calculating actual recombination histories. 



Bound-bound rates 



I'.l 



The spontaneous electric dipole transition rate ^A n ! 
for a nonrelativistic hydrogen atom is given by [79j 

64rr*i£ n , max(U') 



l,l' 

n.n' 



3hc 3 
(i) x U' ( = 



1,1' 

n.n' 



21 + 1 

oo 

x 3 R n >i'(x)R n i(x)dx 



(24) 
(25) 



where e is the charge of an electron, h is the Planck 
constant, and ^X^ n , denotes the radial matrix element 
between the states [n, I] and [n',l'] at order p in the 
multipole expansion. For example, 1 denotes the 
quadrupole rate, and so on. The restriction I' = I ± 1 en- 
forces electric dipole selection rules. Here R n i(x) is the 
radial wave function of an electron in a hydrogen atom, 
with principal quantum number n and angular momen- 
tum quantum number I, at a dimensionless distance x. 
All dimensionless distances are measured in terms of ao. 
For Coulomb wave functions, this integration yields the 
Gordon formula 173: 



(Dv 



LI' 
n.n' 



(-1) 



71 ' —I 



! (n + l)\{n' + l-l)\ 



(Ann' 



(n + n 1 ) 
where V = I — 1 



4 (2Z — 1)! y (n - I - 1)1 {ri - l)\ 
i+i 

7 (n — n ) 



(26) 



/sn+n' — 21 — 2 T , 7 , , n 

W (n,n ,1) . 



W (n, n', I) = 3F1 (u, -n' + 1, 21, w) - 
x 2F1 (v, —n' + 1, 21, w) , 



n + n' 



(27) 



with u = —n + 1 + 1, v = —n + / — 1, and w = 
—Ann'/ (n' — n) 2 . Here 2F1 (a, b, c; x) is Gauss's hyperge- 
ometric function for integer a, b, and c, evaluated using 
the recursion relationship 

(a — c) 2-F1 (a — 1, b, c; x) — a(l — x) [2-F1 (a, b, c; x) 
-2F1 (a + 1, b, c; x)] + (a + bx - c) 2 Fi (a, b, c; x) , (28) 



with initial conditions 



2 F 1 (0,b,c;x) = l, 2 F 1 (-l,b,c;x) = l 



bx 



(29) 
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We use Eqs. (|25p - (|29[) to calculate bound-bound transi- 
tion rates at the beginning of a MLA computation, stor- 
ing them for easy and repeated access. 

We compared the resulting radial matrix elements with 
several values for high n in Ref. (75j and found agree- 
ment to all 3 published digits. We calculated oscillator 
strengths and compared with Ref. (7(| (all transitions 
with n and n' were evaluated, as was the entire Balmer 
series for n < 60) and found agreement to all 6 published 
digits. We also compared with the results in Ref. [77} (in 
which oscillator strengths were computed up to n = 500 
for An < 5) and found agreement to 5 digits. We at- 
tribute the difference in oscillator strengths to the fact 
that a polynomial expansion of 2F1 was used in Ref. [77} , 
rather than the more stable recursion relationship. We 
also compared with the dipole one-photon rates used for 
the n max = 30 MLA computation of Ref. [5l|. Most 
rates agreed to 7 or more significant figures. Transition 
rates between s and p orbitals only agreed to ~ 5 signif- 
icant figures. We ran our MLA model using the rates of 
Ref. [51J and verified that these small disagreements do 
not lead to any differences in x e (z) at the desired level 
of accuracy. Given the high quantum numbers consid- 
ered, it was important to verify that no numerical in- 
stability plagues our numerical implementation of these 
recursions. We thus checked matrix elements computed 
using Eqs. (l2"T5)) - (j2"gj) against values estimated using the 
WKB approximation, as detailed in the Appendix. 



The initial conditions of the recursion are [75 



G n,o =\ln7KZ 7T7 ( 4n ) e 



n -2n 



G 



2 (2n-l)! 

^ g2n— 2k~ 1 atari (tik) 



Vl-e-^ (l + n 2 K 2 ) 
G™" 2 '"- 1 = (2n - 1) (1 + nV)^" 1 '" 
1 + n 2 K 2 



n+2 w n,0 



/-in 



G 



n — 1 , n — 2 



In 1 



(32) 



These matrix elements are tabulated at the begin- 
ning of each MLA run for all I < n < n max , and 
10~ 25 < n 2 n 2 < 4.96 x 10 8 ; this range of k is parti- 
tioned into 50 logarithmically spaced bins, with each bin 
containing 11 equally spaced K values. Bound- free matrix 
elements were compared with tabulated values for low n 
in Ref. [78[ and agreed to all 4 listed digits. Matrix ele- 
ments were also compared with those used in Ref. [5lj |; 
we found agreement to one part in 10 7 , aside from s — p 
transitions, as already discussed. 

The recombination rate to [n, I] as a function of energy 
is then 

«m (E e ) = 4 /™ Yl max {I, I'} 9%, (33) 



3n 2 (A;Tm) 



3/2 L^i 

i'=i±i 



with 



2. Bound-free rates 



Ql l' = I I 



n 2 E, 



U' 

■V 1 1, . K 



(34) 



Bound-free rates are evaluated using the same princi- 
ple, but one of the two states used to evaluate matrix 
elements must be a continuum Coulomb wave function. 
The resulting matrix element is [801 ] 



9n,' K =-2l x i R n i{x)F K [ / (x)dx, 



1 



(30) 



where F K i' is the continuum Coulomb wave function 
for a recombining photoelectron with angular momen- 
tum quantum number I' and dimensionless energy k 2 — 
Ee/I-R = 7^ — 1/n 2 . The energy of the outgoing photon 
is hv. This integral may also be evaluated in terms of hy- 
pergeometric functions, which in turn yields a recursion 
relationship for g\^ K [781 : 



G 



1,1' 



1,1' 



(2^VCTin:=o(i+A 2 ) 

G 1 ^ 1 - 1 = [4 (n 2 -l 2 )+l (21 - 1) (1 + nV)] 



(n+l)\ 



An 2 (n 2 -l 2 ) 1 + (I + 1)V 



G 



1,1+1 



G,, 



= [4(n 2 -l 2 )+l(2l + l)(l + n 2 K 2 )]G^ 1 
- An 2 \n 2 -(1 + if] (1 + l 2 n 2 ) Gl+t l - (31) 



At each value of Tm, the tabulated matrix ele- 
ments, Eqs. ((6|) and ([33| are used to calculate thermally 
averaged recombination rates, using an 11-point Newton- 
Cotes 81] formula for the integration and neglecting 
stimulated emission. Large bins are added into the 
integral until it has converged to a fractional precision 
of 5 x 10~ 15 . We com par ed our values with integrated 
rates tabulated in Ref. [78| and found agreement to all 4 
listed digits. Comparing with the rates used in Ref. [5lj ]. 
we found agreement to one part in 10 7 , aside from s-p 
transitions. 



In Saha equilibrium, 

rila nl (E e ) [l+N(E 7 ,T R )} P M (E e ,T M ) 

= m x n,lN (E-f,TR) /3 n i (E 7 ) , 

and so by the principle of detailed balance, 

dE P J nl (£ 7 ) = ME [dE e a nl (E e ) 

Xn,l J 

[l+M(E 7 ,T n )] 



(35) 



AA(^ 7 ,T R ) 



Pm(E £ 



(36) 



We verified that our computed thermally averaged re- 
combination and ionization rates satisfied this equality 
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to machine precision. We also checked bound-free ma- 
trix elements computed using Eq. (|31l) against values es- 
timated using the WKB approximation, as detailed in 
the Appendix IVI1 



B. Sparse-matrix technique 

The key to making the recombination problem 
tractable for high values of n max is the sparsity of 
Eqs. (TlH) and (IT51) . Dipole selection rules only allow cou- 
pling of states with angular momentum quantum num- 
bers I and I' if V = I ± 1. It is easiest to understand how 
sparsity simplifies the problem with a slight change of no- 
tation. We can compose the vector x (with components 
x n ^i) of excited-state populations, as 



(37) 



where l max = n max — 1 and vi denotes a vector of the pop- 
ulations of all states with angular momentum /, except 
for the Is state. Specifically, 



Vj 



where 



Sn min +l,i 



2 if I = 0, 
l + l if I ^ 0. 



(38) 



(39) 



The source vector s can similarly be written by con- 
catenating source vectors s/; each s; feeds all states with 
angular momentum I. 

The rate matrix may be similarly built of submatrices 
M/,c , as illustrated in Fig. [T] The complete rate matrix 
is block tridiagonal, and the blocks decrease in dimension 
as I increases. The matrix Mj ;/ has components 



Kf = C Uu+™+ E <>>n> ~ T n'n>- (40) 



In the steady-state approximation, Eq. (|14[) can be 
rewritten as a system of matrix equations. If Z = 0, 

Mo,otfo+Mo,i3i = 3b- (41) 

If < I < z max , 

M M _i#i_i + Mtjvt + M M+1 ?7 /+ i = si. (42) 

To close the system, we must truncate the hierarchy 
by excluding states with n > n max as both sources 



and sinks, which is equivalent to setting A 
max{n, n'} > ri max . Then for I = l max , 



i±i 



M, 



-liiL-i + Mi. 



for 



(43) 



It might be possible to approximate the correction due 
to this truncation error, using asymptotic expressions for 



and Saha equilibrium abundances for n > n B 



This will only work if n max is sufficiently high for nearly 
perfect equilibrium Saha equilibrium to hold between 
states with n > n max and the continuum. 

At any given time step, the actual quantity of interest 
is not the inverse T _1 of the rate matrix but the solution 
set {vi} to the steady-state rate equations. The closed 
form solution to Eqs. l|iT |) -([i3 | is 



vi 



K, 



si 



M u+1 v l+1 + 



^max, then 



1-1 

l'=0 



S/,/'S// 



(44) 



Here 

K ; = 
and 



M 



K, 



si 



E(-i) r " 

i'=0 



Sli'Sv 



00 



M z ,j_iKi_iM. 



i-i,i) 



if I = 0, 
if I > 0, 



Mj, w K w if i = 1-1, 
S u+1 M 1+M K t if i < I- 1. 



(45) 



(46) 



(47) 



Our new MLA code, RecSparse, operationally imple- 
ments this solution at each time step as follows: 

1. Using the values of Tr and x e , Tm is calculated 
using the results of Sec. Ill Bl 

2. All relevant Mi^/ and s*/ are computed using 
Eqs. (gDD and (TTH)) and stored. 

3. All K ; and S M are computed using Eqs. (BH]) - (fT7|) 
and stored. 



4. Equation (I45|) is applied to obtain the solution for 
i?U ■ 



5. Equation 
all vi. 



is iterated to obtain the solutions for 



The free-electron fraction x e is then evolved forward in 
time using {vi } and Eq. ([20)) . It would also be interesting 
to compute the cumulative spectral distortion emitted by 
the line and continuum processes responsible for recom- 
bination [43|, [Ml Is2] ■ This fractional perturbation of 
1CP 7 to the blackbody intensity of the CMB could be de- 
tectable with future experiments and would offer a test 
both of our understanding of recombination and of new 
physics behind the surface of last scattering (e.g., time 
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FIG. 1: Schematic of the sparse rate matrix T with components given by Eq. l|15f) and submatrix building blocks given 
by Eq. ()40f) . Boldface zeroes denote block matrices of all zeros, and enforce the dipole selection rule that the initial state 
I' angular momentum obeys I' = I ± 1, where / is the final state angular momentum. The submatrix M ;i / has dimension 
(imax — ^min + 1) x (ri max — n' mln + 1), where n m i n = 2 if I = 0, and n m i n = I + 1 if I > 0. Note that submatrices Mij on 
the block diagonal of the larger rate matrix T are themselves diagonal, as seen from Eq. (|40p and the fact that in the purely 



radiative case, F 



if n ^ n' and I = I'. 



variation of fundamental constants, energy injection by 
decaying/annihilating dark matter) [8314871 ] . This and 
the development of a fast code for Planck data analy- 
sis including all the relevant physical effects will be the 
subject of future work. 



Numerical methods 



RecSparse begins at z = 1606, assuming Saha equi- 
librium to compute the initial value of x e and setting 
Tm as discussed in Sec. Ill Bl Excited-state populations 
are obtained using the method of Sec. IIIIBI Subma- 
trix inversions are implemented using the double preci- 
sion routine DGESVX from the LAPACK library [H]. 
Time evolution of x e (z) with Eq. ([2T))) is implemented us- 
ing the 5 th -order Runge-Kutta-Cash-Karp (RKCK) im- 
plementation in Numerical Recipes [891 ]. The rapid time 
scale for return to Saha equilibrium introduces a stiff 
mode into the equations at early times, necessitating care 
in the choice of a stepsize for the integrator. We were 
able to achieve relative precision of e ~ 10~ 8 by placing 
59 time steps at z > 1538 and 250 steps in the range 
200 < z < 1538, partitioning each interval into equally 
sized steps in A In a; relative errors were estimated by 
halving step size and comparing values of x e (z) at iden- 
tical time steps. The computation time t comp for Rec- 
Sparse scales as t CO mp oc n^, where 2 < a < 3. This is 
an empirical estimate for the range of n max that we have 
explored, and may not extend to higher n max values. In 
contrast, for standard MLA codes, t comp oc n® lax . We 
can calculate recombination histories for n max = 200 in 
4 days on a standard workstation; this would likely take 
weeks using a conventional MLA code. 



IV. EXTENSION TO ELECTRIC 
QUADRUPOLE TRANSITIONS 

Early work on recombination highlighted the impor- 
tance of forbidden transitions, as half of the hydrogen 
atoms in the Universe form by way of the 2s — > Is "for- 
bidden" transition (i^, |4l| . Recent work has included ad- 
ditional "forbidden" transitions in the MLA treatment, 
namely, two-photon transitions (ns —¥ Is and nd — > Is) 
in H [48-51], two-photon and spin-forbidden transitions 
in He 52H55 J. a s well as electric quadrupole (E2) transi- 
tions in He|^,|73|. 

Until this work, the impact of E2 transitions in H on 
recombination has not been considered, even though they 
are optically thick for transitions to/from the ground 
state. For optically thick lines, the overall transition 
rate is proportional to A^/t^,. Since t^[, cx A 1 ^,, 
the overall transition rate is independent of the rate co- 
efficient. Transitions such as electric quadrupoles, which 
seem "weaker" judging from rate coefficients alone, can 
thus be as important as "stronger" transitions, like the 
Lyn lines. For example, this is why the semiforbidden 
He I 59lA line is important in cosmological recombina- 
tion [69|, [70j • We thus include E2 quadrupole transitions 
in our MLA computation to properly assess their rel- 
evance for cosmological recombination. Ml (magnetic 
dipole) transition rates in H are typically suppressed by 
an additional factor of 10 — 10 8 , and are thus negligible 
90]. 



A. Rates 

The electric quadrupole (E2) Einstein A-coefficient for 
transitions from states [n,l] to states [n',1 1 ] is [91| : 



(2) A l',i 



I5g a c 4 



<nZ||Q( 2 )||n'0 



(48) 
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where the quadrupole matrix element is 

(nl\QM\n'l') = {l\\CW\\l') ^X%. (49) 

The matrix elements of the reduced angular tensor op- 
erator are given by 



(i\\c i2) \\n = (-i) l VWgf 



I 2 V 

ooo 



(50) 



where the last factor is the well-known Wigner-3 J sym- 
bol. This operator is defined as 



(l\\C^\\l') = (-1 



/ ,n I I k V 

-771 q w! 



x ^^ TJ (lm\Y kq (e, ( f > )\l'm').(51) 
The dimensionless radial quadrupole integral is 

/•OO 

(2) *n'*» = / x i R n , v {x)R n i{x)dx. (52) 



The radial matrix element for the nd — > Is transition is 
a special case of Eq. (B.13) of Ref. 92] with 77' = 1: 



WX°* = (-ir- x 2 6 n 4 



(77 + 2)! 
(n-3)! 



1/2 



( n - u 

I i -i \n+4 

(n + 1) 



(53) 



the transition nd — > is is highly probable to be immedi- 
ately followed by a transition Is — > np. This yields a net 
77tf — > np transition, analogous to an /-changing collision, 
which occurs with forward rate = x nd ^A { 2 n . The 

reverse process occurs with rate ^r"'^ = x np ^A^ n D, 
where D is a factor relating forward and backward rates. 
If the p and d states were in equilibrium, the two rates 
would cancel, so by the principle of detailed balance, 
D = (x n d/x np ) c(i = 5/3, where "eq" denotes an equi- 
librium value. The net np nd transition rate due to 
E2 transitions is thus 



Xnd 



0.2 



r 



n p 



(56) 



Since this overall rate obeys the Al = ±1 selection rule, it 
can be numerically implemented within the same frame- 
work as the dipole rates. 



RESULTS 



We ran the RecSparse code for a variety of n max val- 
ues. Here we omitted E2 transitions to focus on the effect 
of deviations from statistical equilibrium and increasing 
n mn . We begin by discussing deviations from equilib- 
rium, and proceed to discuss the recombination history 
and numerical convergence with n max . 



B. Inclusion in multilevel atom code 

The obvious way to include quadrupole transitions into 
our MLA code would be to generalize Eq. (l4"2"j) to include 
Al = ±2 transitions: 



M M+2 v /+2 + M M+ iv i+ i + M M v ; 

+ M u _iV;_i + M U _ 2 V/_2 = Sj. 



(54) 



The resulting system is obviously not as sparse as in the 
dipole case, and solving for all v/ would be computation- 
ally more expensive, slowing down the whole MLA code. 
Since the contribution from even the largest quadrupole 
rates may turn out to be small, we pursue a computa- 
tionally less expensive approach. 

Higher energy E2 transitions will proceed much faster 
than lower energy ones, since E2 rates scale as u!^ n ,. In 
particular, transitions to and from the Is ground state 
will dominate any other quadrupole contributions to the 
recombination problem, since 



(2) A - 2 



, .5 



q 2 (n 2 



1 



;> 10 3 if q > 2. (55) 



Moreover, the nd — > Is lines are optically thick for small 
n. We thus restrict our consideration to nd o Is transi- 
tions, since other quadrupole transitions are "corrections 
to a correction." A further simplification follows if we re- 
call that the Lyn lines are all optically thick [5l|. Thus, 



A. State of the gas 

The assumptions of statistical equilibrium between dif- 
ferent I sublevels within the same n shell and Boltz- 
mann equilibrium between different n states fail at late 
times, as discussed in Sec. IIIII Furthermore, as reac- 
tions become inefficient on the Hubble time scale and 
x e (z) freezes out, Saha equilibrium between the contin- 
uum and excited states of H also fails. Below, we discuss 
each of these failures quantitatively. 



1. Populations of angular momentum sublevels 

At early times, the populations of hydrogen atoms in 
states with the same n but different angular momentum 
I are in statistical equilibrium [see Eq. ([HJ) ]- Radiative 
transitions do not include reactions that are I changing 
but 77 conserving. The I sublevels must thus be kept in 
equilibrium by a combination of sequences of allowed ra- 
diative transitions and atomic collisions. These processes 
become inefficient at later times, leading the different I 
sublevels to fall out of equilibrium. Both the TLA treat- 
ment of Peebles and the later MLA treatment of Sea- 
ger et al. rely on the statistical equilibrium assumption 
IS EH- Our RecSparse code relaxes this assumption 
and follows the populations of all I sublevels separately. 
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For n > 5, the resulting populations are marked by 
several features, shown in Figs. [2] and [3] at early and late 
times, respectively. We use 

Aln,! = XnJ - X° q ; (57) 

to compare actual with equilibrium populations, where 
Co (2Z + 1) 



(58) 



< 



Deviations begin modestly at early times ( Ax n< i/x 

0. 1% for 1300 ^, z <, 1600) but are quite large by late 
times ( Ax n .i/x c ^i ~ 60% by z <, 600). 

Lower I states depopulate efficiently, and are signifi- 
cantly underpopulated relative to statistical equilibrium 
expectations. States with I = can only make downward 
dipole transitions in n if I' = 1. These rates are several 
order of magnitude lower than Lyman-series rates with 
the same An, and so I — states depopulate less effi- 
ciently than other low-/ states. This explains the upturn 
at the lowest I values. The Al = ±1 selection rule im- 
plies that higher I states couple efficiently to neighboring 
bound states (V = I ± 1) with a limited range of acces- 
sible n', since n' > I'. These states thus depopulate less 
efficiently than states with lower I due to this bottleneck. 

The recombination rate a n i peaks in the range 0.3 ^ 
l/lmax & 0.4. Together, these facts imply the presence 
of a peak in As„ : j/j;^ q ; , which turns out to occur in the 
range 32 <5 I <, 37 for a wide range of n at all times. The 
transition to x nt i/x c ^ l > 1 occurs in the range 16 ^ I ^5 
21, also for a wide range of n at all times. At very high 

1, recombination rates are so slow that these states are 
again underpopulated relative to statistical equilibrium, 
though less dramatically than they are at low I. 

The observed amplitude and shape of the curves in 
Figs. EE ] qualitatively agree with the results in Refs. 
60]-|61], including the upturn near the lowest I and sharp 
minimum at I = 2. The minimum is due to fast Balmer 
transitions out of the I = 2 state. When we computed a 
recombination history with these rates {nd — > 2p for n > 
2) artificially set to zero, the minimum moved to / = 1, as 
shown in Fig. [4] It is interesting that the curves in Figs. 
[2H3] exhibit the same behavior with I as the departure 
coefficients of Ref. j93[ , which describe neutral hydrogen 
(also in the steady-state approximation) in interstellar 
H II regions. 

RecSparse only takes into account radiative transi- 
tions, and omits / and n-changing collisions. These rates 
would flatten all the curves in Figs. lessening de- 

viations from statistical equilibrium between the differ- 
ent I sublevels [6l|. Indeed, the assumption of statisti- 
cal equilibrium between these states at all times is for- 
mally equivalent to the limit of infinite Z-changing colli- 
sion rates. Theoretical estimates for collisional rates all 
depend on different assumptions and tabulated rates dis- 
agree by factors of two or more (see, e.g., Ref. [H]]). As 
a function of redshift z, we estimate the ratio /£° n = 



(in %) 



-0.02 



-0.04 



-0.06 




1+1 



FIG. 2: Early time deviations from statistical equilibrium be- 
tween different I at fixed n and n max , as computed by Rec- 
Sparse. 



ijaXe(lnitni of collisional to radiative transition rates out 
of the state [n,l], where q n i is the collisional rate coeffi- 
cient (in cm 3 s _1 ), and t n \ is the total radiative lifetime of 
the state, including stimulated emission and absorption. 

Using the rate coefficients in Ref. J95|, we estimate 
that collisional rates (per unit time) are of the same order 
of magnitude as radiative rates for n <; 52 at z ~ 1600, 
n ;> 83 at z ~ 1080, n <; 160 at z - 740, and n ^ 250 
at z ~ 200. In other words, as the primordial gas cools, 
collisions come to only influence the highest H energy 
levels, which contain the least bound electrons. This 
occurs because of the exponential decrease in the free- 
electron density r]^x e in the early stages of recombina- 
tion, which drives down collision rates accordingly. Near 
z ~ 1600 and shortly thereafter, radiative rates alone are 
high enough to keep the excited states in Z-equilibrium. 
Collisions thus have little effect on x e (z) at early times. 
There may, however, be a window at some intermediate 
redshift, when collision rates are still relatively high, but 
departures from Z-equilibrium are large enough to war- 
rant including collisions in the recombination model. A 
full calculation is necessary to understand the actual ef- 
fect. A final answer on the effect of resolving I sublevels 
on both the recombination history x e (z) and the recom- 
bination spectrum awaits a robust theoretical calculation 
of the relevant collisional rates. This is an area of future 
investigation. 



2. Populations of Rydberg energy levels 

We may also compare the total population of the n th 
energy level to values in Boltzmann equilibrium with n = 
2: 
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FIG. 3: Deviations from statistical equilibrium between different I at fixed n and n max , shown as computed by RecSparse at 
a variety of times through the recombination process. The left panel shows results for states with n = 25, while the right panel 
shows results for states with n = 140. 
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at fixed n and n max are shown at a variety of times through the recombination process. The left panel shows standard results 
with RecSparse. The right panel shows the results obtained if I — 2 Balmer rates are artificially set to zero in the code. This 
figure highlights the relatively rapid 1 — 2 Balmer transitions as the origin of the / = 2 dip. 



The recombination rate to states with n > 2 is greater 
than the downward cascade rate, creating a bottleneck 
to depopulating these states. This bottleneck causes an 
over-population of the excited states compared to the 
equilibrium values of Eq. ([5^| . as shown in Fig. [SJ The 
ratio x n /Xn° Uz is O (1) at early times but grows as high 
as 3 x 10 4 by z = 555. The ratio approaches a constant at 
high n, as energy levels get closer to the continuum and 
the energy differences between successive levels shrink. 

Relative to n = 2, excited states are over-populated, 
but there is no population inversion or cosmic maser. Ex- 



cited states are still less populated than the n = 2 energy 
level, just not as dramatically as they would be if Eq. (|59|) 
held. Among highly excited states, some pairs of levels 
do exhibit population inversion. For effective maser ac- 
tion, inversion must occur between pairs of radiatively 
connected levels, and the coherence of the radiation field 
must not be destroyed along the line of sight. This effect 
will be explored in detail in future work. In extremely 
dense structure-forming regions, more dramatic popula- 
tion inversion may result and lead to local masing; if 
these masers were observed, they could offer interesting 
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new probes of structure formation near z ~ 1000 as well 
as the physics of reionization [96|. 

Recombination becomes inefficient at late times; i.e., 
the recombination time scale [a-B{T)x e n-£{\~ 1 becomes 
longer than the age of the Universe. Saha equilibrium 
expressions for x e and X\ s fail dramatically at late times. 
The free-electron fraction x e freezes out and is higher 
than the Saha equilibrium value, and thus x\ a is lower 
than the Saha equilibrium value. Excited states are 
overpopulatcd relative to the ground state, but still not 
enough to be in Saha equilibrium with the continuum. 
The tower of excited states is thus also underpopulated 
relative to Saha equilibrium, as shown in Figs. [6] and 
[7J Lower energy levels fall out of Saha equilibrium faster 
than higher energy levels. Higher energy levels are clos- 
est to Saha equilibrium, but at late times (z ~ 200), 
even the population of the n — 250 level is nearly 10% 
below its Saha equilibrium value. Modeling the effect 
of states with n > n mayi may require Saha equilibrium 
abundances to hold in the regime past the cutoff. To this 
end, it is important to properly model atomic collisions 
(which would push atoms towards Saha equilibrium at a 
lower transitional value of n max ), and apply even greater 
computational resources to obtain x e (z) for even higher 
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FIG. 7: Actual population of energy shells compared to Saha 
equilibrium values, shown for several n values as an explicit 
function of cosmological redshift z. 



n max = 4, 8, 16, 32, 64, 128, and 250. We define a relative 
error: 



Axl (z) = £e max (z) 



(60) 



B. The effect of extremely high-n states on 
recombination histories and the CMB 



To explore the relative convergence of x e (z) over a wide 
logarithmic range of n max values, we computed x e (z) for 



Here n l max is the i th n max value. We show the result- 
ing recombination histories and Ax z e (z) in Fig. [H As 
"max increases, the larger number of pathways to the 
ground state makes recombination more efficient, de- 
creasing Xe max (z) and making Ax e (z) positive. The rel- 
ative error Ax l e (z) shrinks with n ma x, indicating that 
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relative convergence is taking place, as demonstrated in 
Fig. [§] Note, however, that the relative error may not 
be a good proxy for the absolute error. Suppose that 

the absolute error is given by i" m " = Axf 38 ' 1 + x e , 
where Ax| bs ' 1 = A(ro max ) P , for some normalization A 
and power-law index p < 0. Then it is easy to show 
that for n max = 2t&£, Axl/Axf s ' ! = (1 - 2*>). In 
other words, the relative error will underestimate the ab- 
solute error. To demonstrate absolute convergence, one 
should demonstrate that the physics neglected by ignor- 
ing transitions to n > n max docs not cause large changes 
in x e (z). We also calculated recombination histories for 
n max = 20, 50, 90, 105, and 160. 

We may also assess the effect of the computed changes 
in x e (z) on the CMB Cg's. To this end, we replace the 
usual table generated and used in the RecFast module 
of CMBFast with a table of our own output for dif- 
ferent n max values, smoothly stitching our history onto 
the usual RecFast history at the boundaries z = 1606 
and z = 200. We tried a variety of smoothing schemes 
including no smoothing at all, and determined that the 
resulting error was at most 10% the change already in- 
duced by varying n max . The choice of smoothing scheme 
is thus a "correction to a correction" and does not al- 
ter the conclusions of our analysis. In particular, the 
number of sigmas at which power spectra corrected and 
uncorrected for higher-n levels can be distinguished will 
change by at most 10% of itself as a result of changing the 
smoothing scheme. The statistical significance of higher- 
n shells will thus be essentially unchanged by the choice 
of smoothing scheme. The results for temperature and 
E-mode polarization anisotropy power spectra (Cj T and 
Cf E ) are shown in Figs. ITUl and fTTl respectively. Here 
we also define a relative error: 



AC**'* = cf x '< 



-C. 



xx, < 

i 



(61) 



Here XX denotes the TT or EE label of the power spec- 
trum under consideration. The relative error AC XX:1 is 
always positive, indicating that increasing rt max also in- 
creases Cf x , as shown in Figs. 111)1 and ITT1 The common 
(TT and EE) origin for this effect is clear from Fig. [5] 
Higher n max makes recombination more efficient, driving 
down the freeze-out value of x e (z) and the residual opti- 
cal depth r, leading to the high-/ plateaus seen in Fig. [JO] 
and [TTJ As a result, the smearing out of primary CMB 
anisotropies by relic free electrons, Ci — > Cie~ 2T [35j, is 
less dramatic when rt max is increased. The relative error 
AC £ ' l shrinks with increasing n max . 

Taken as a proxy for the absolute error, AC XX '' may 
be compared to a crude (cosmic variance) estimate of the 
required accuracy of C xx predictions in the damping tail: 



AC, 



xx 



C 



XX 



3 x 10~ 4 f , 1/2 



'sky 



(62) 



Here / s k y is the fraction of the sky covered by a CMB 
experiment. For / s k y = 0.70, results are shown in Figs. 



[TOl and [TT] and we see that only for n max = 250 does 
the relative error shrink to a level comparable with the 
cosmic variance. The ultimate aim is for the total correc- 
tion from recombination physics to be less than statistical 
errors, so any individual contribution such as the trunca- 
tion error at n max should be -C la. In any case, collisions 
must be properly included to show absolute convergence, 
and so this should be a key focus of future work on highly 
excited states in hydrogen recombination. To more real- 
istically assess the importance of high-n states, ACf x 
should be compared with a realistic error estimate for 
Planck. 



C. Statistical significance of corrections to 
recombination history 

As a test of the importance of the modified recom- 
bination history for Planck, we have compared our cor- 
rections to the power spectrum ACi with the forecast 
Planck error bars. The comparison is done by means of 
the statistic 



Z= J^FwACiACv 



(63) 



where Fu/ is the Fisher matrix for the CMB power spec- 
trum. For the temperature-only case, £ ranges from 
2 to £ max and hence F is an (£ max - 1) x (£ max - 1) 
matrix; when polarization is included, F expands to a 
3(^ ma x — 1) x 3(£ max — 1) matrix incorporating TT, EE, 
and TE spectra. The Z statistic is the number of sig- 
mas at which the corrected and uncorrected power spec- 
tra could be distinguished assuming perfect knowledge 
of the cosmological parameters, and hence represents the 
largest possible bias (in sigmas) on any combination of 
cosmological parameters in any fit that incorporates the 
CMB [5JJ. We use the forecast noise and beam curves 
for Planck data 70 GHz (Low- Frequency Instrument) and 
100 and 143 GHz (High-Frcqucncy Instrument) channels 
in the Blue Book [22|, and assume a usable sky fraction 
Of /sky = 0.7. 

The computation considering the difference between 
the n max = 128 and 250 curves gives a Z value of 0.36. 
However, the actual error in the n max = 128 calculation 
is somewhat greater because even the n max = 250 cal- 
culation is not completely converged. If the error in the 
CiS scales as ~ njj,^ and has a shape that varies slowly 
with n max , then our value of Z should be increased by a 
factor of [1 - (250/128) p ]~ 1 ; for p « -1.9 (as suggested 
by Fig. \§§ this is 1.39. Thus if the power-law extrapo- 
lation is to be trusted there is a 0.50<r error (Z = 0.50) 
in the CMB power spectrum if one restricts attention to 
"■max = 128, and a ~ 4 times smaller error [Z = 0.14) at 
"■max = 250. A similar comparison between n max = 64 
and 250 implies an error of Z = 1.79 at n max = 64. This 
suggests that in the purely radiative problem the CMB 
power spectrum is converged (in the sense that our re- 
maining errors are small compared to projected Planck 
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errors) at n max > 128; however this issue will have to be 
reconsidered in future work when collisions are included. 



D. The effect of electric quadrupole transitions on 
recombination histories and the CMB 

Using the treatment of Sec. HVI and an integration step- 
size fine enough to obtain a fractional accuracy of 10~ 10 
in x e , we compute the effect of E2 quadrupole transitions 
on cosmological hydrogen recombination for several val- 
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FIG. 10: Relative errors between temperature anisotropy 
spectra Cj T computed using CMBFast, modified to include 
successively more accurate RecSparse recombination histo- 
ries. Pairs of n max values used for the comparison are indi- 
cated in the legend. Cj T increases with n max , as discussed 
in Sec. IV Bl The correction shrinks with increasing n max . 
The long dashed line indicates the cosmic variance target for 
ACi/Ce, as discussed in the text. 



ues of n m ax- We can parametrize this effect using 



Ax e 



■ E2 transitions e lwith E2 transitions 



(64) 



and 



ACe — C£| with E2 transitions ^Ino E2 transitions ■ (65) 
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FIG. 11: Relative errors between E-mode polarization 
anisotropy spectra Cf E computed using CMBFast, modified 
to include successively more accurate RecSparse recombina- 
tion histories. Pairs of n max values used for the comparison 
are indicated in the legend of Fig. 1101 Cf E increases with 
n-max, as discussed in Sec. IV Bl The correction shrinks with 
increasing n max . The long dashed line indicates the cosmic 
variance target for ACe/Ce, as discussed in the text. 



Note that unlike the case of varying n max , these are the 
absolute errors induced by ignoring E2 transitions. 

The results are shown in Fig. [T^l The maximum effect 
of E2 transitions occurs at z ~ 800 with a fractional en- 
hancement of Ax e /x e ~ 10 , and the calculation seems 
well converged by n max = 30. Corrections due to higher 
excited states would be a correction to a correction, and 
so we ignore them. Although the correction from E2 
transitions is small, it extends over a broad epoch at late 
times after reaching its maximum. To determine if this 
could affect CMB anisotropics in an observable way, we 
modify and run CMBFast [97J using recombination histo- 
ries computed with/without E2 transitions. We incorpo- 
rated RecSparse recombination histories including E2 
transitions into CMBFast by applying the same method 
employed in Sec. IV Bl 

Running the recombination histories including E2 
quadrupole transitions through CMBFast gives a max- 
imum change ACe/Ce ~ 3 x 10~ 6 in both temperature 
and polarization, negligible compared to cosmic variance. 
Thus E2 transitions in hydrogen are negligible for CMB 
applications. 



VI. CONCLUSIONS 

We have developed a new recombination code, Rec- 
Sparse, optimized for tracking the populations of many 



FIG. 12: Fractional difference between recombination histo- 
ries with/without E2 quadrupole transitions included for dif- 
ferent values of n ma x- The net effect is always to speed up 
recombination. 



energy shells of the hydrogen atom while resolving angu- 
lar momentum sublevels. The code runs more quickly 
than would be anticipated using simple scaling argu- 
ments, which would yield the the scaling t comp cx n^ ax . 
Using RecSparse, we find empirically that for the range 
of n max values used, computation time scales as t comp cx 
n maxi where 2 < a < 3. With this code, we have com- 
puted cosmological hydrogen recombination histories for 
a series of n max values going as high as n max = 250 and 
explored the highly nonequilibrium state of the resulting 
atomic hydrogen. 

The resulting correction Ax e (z) satisfies 
Ax e (z)/x e (z) < 0.01 for z > 200 when n max = 250 and 
converges with Ax e (z)/x e (z) cx The correction 

to the C/s becomes of order the cosmic variance 
when n max = 250. In light of realistic error estimates 
for Planck, the resulting CMB anisotropy spectra 
Cf x are converged to 0.5c at Fisher-matrix level for 
«ma X = 128 in the purely radiative case, assuming error 
extrapolations may be trusted. 

To definitively answer the question of absolute con- 
vergence, collisions must be included to speed the ap- 
proach to Saha equilibrium at high n, allowing a conclu- 
sive treatment of states beyond the truncation limit, with 
n > n max . Future work should also properly account for 
the overlap of the Lyman resonance line series at high 
n. It will also be interesting to determine if there is co- 
herent stimulated emission between excited states, given 
its relevance for the detectability of faint CMB spectral 
distortions from the epoch of recombination. Finally, the 
sparse-matrix methods applied here or similar techniques 
could be profitably applied in the development of fast re- 
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combination codes for CMB data analysis, even at early 
times in recombination, when only lower values of n max 
are relevant. 
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Appendix: WKB approximation for radial dipole 
integrals 

The development of laser spectroscopy of high-n states 
in hydrogen and other atoms, along with the study of 
nonlinear and multiphoton ionization, required the com- 
putation of dipole radial matrix elements for high and 
even fractional quantum numbers in a Coulomb or per- 
turbed Coulomb potential 98]. Until adequate algo- 
rithms for these computations were ultimately developed, 
the Wentzel, Kramers, Brillouin, and Jeffreys (WKBJ) 
semiclassical approximation (quite accurate for n 1) 

[H HI proved a useful tool for estimating (n) A^V„. At 
high n, radial wave functions in the Coulomb potential 
have a large number of nodes and thus a short wave- 
length A. For the WKB approximation to be valid, it is 
necessary that \d\/dx\ <C 2ir. Because of the large num- 
ber of nodes in the Coulomb wave functions at high n, 
the WKB approximation is ideally suited to estimating 
matrix elements for transitions between high n. 

In the classically allowed region, the nonrelativistic 
WKB radial wave function for a hydrogen atom is 



with 
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where the inner classical turning point x\ is a solution 
of the equation k n i(x) — 0. Substituting Eq. (IA2[) into 
Eq. (|25[) for the dipole matrix element, and making sev- 
eral additional approximations, the following expression 
is obtained if \n' — n\ <C n, n' and n, n' 3> I [99J : 
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with s = n - n', Al = V - I, l c = + + 1) /2, 
n c = Inn' I (n + n'), and e 2 = 1 — (/^/n 2 .) ■ Here e is 
the eccentricity of a Keplerian orbit with the quantum 
numbers n c and l c , and J s (x) is a Bessel function of the 
first kind. These estimates agree with matrix elements 
computed using Eq. (|2"BT) to a precision of 5%-50%; the 
agreement worsens as \n' — n\ — > n, n' 



If / <C n' , n and 



then pot 



72 

(l) v l±l,l o 1 ( 2 _x 

1 'X„, „ = 2 y= (nn) y 



A" 



2/3 



6 



T^l/3 



l^y 
6 



(A4) 



with y = |ri~ 2 — n' 2 \. Here K s (x) is a modified Bessel 
function of the second kind. These estimates agree with 
matrix elements computed using Eq. (|26[) to a precision 
of l%-20%; the agreement worsens as s shrinks, at which 
point Eq. (|A3j) becomes more accurate. 

A WKB estimate of bound-free matrix elements is ob- 
tained by making the substitution n' i/K in Eq. (|A4|) 
[lOOj . The resulting estimate is reasonable if I <C n^r 1 
and agrees with matrix elements computed using Eq. 
(f3"Tj) to a precision of 50%. This analysis confirms that 
the high n and I values under consideration do not afflict 
our evaluation of Eqs. (|2l?|) or (|3~Tj) with any instability 
that would throw computed rates off by orders of mag- 
nitude. 
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